Lenosky's energy and the phonon dispersion of graphene 
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We calculate the phonon spectrum for a graphene sheet resulting from the model proposed by T. 
Lenosky et al. (Nature 355, 333 (1992)) for the free energy of the lattice. This model takes into ac- 
count not only the usual bond bending and stretching terms, but captures the possible misalignment 
of the p z orbitals. We compare our results with previous models used in the literature and with 
available experimental data. We show that while this model provides an excellent description of the 
flexural modes in graphene, an extra term in the energy is needed for it to be able to reproduce the 
full phonon dispersion correctly beyond the F point. 

PACS numbers: 81.05.Uw, 63.22.-m, 63.20.D 
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I. INTRODUCTION 

The phonon dispersion of graphite is a reccurrent topic 
in the literature, and the last years have seen a revival 
of interest due to the successful isolation of a single layer 
of graphite, graphene 1 . Due to the very weak interaction 
between the planes, the phonon spectra for graphene and 
graphite are essentially the same, except at very low fre- 
quencies where the splitting of the out of plane acoustic 
mode in graphite is noticeable^. In general the theoretical 
models most successful in describing phonon dispersions 
are first principle ones, and graphite is not the exception. 
However simple, analytical models that would give a good 
qualitative description of the system are highly desirable. 
These models are valence force models and usually re- 
quire many free parameters to give accurate results. The 
generally accepted as best fitting valence force model for 
graphite is the one given in Ref. [3j, in which one con- 
siders interactions up to four nearest neighbors and 20 
fitting parameters. More recently a five nearest-neighbor 
model was shown to give a very good fit to new available 
experimental data-. Also a valence force model based 
in the symmetries of the lattice was proposed^. Simpler 
valence force models in which only nearest neighbor in- 
teractions are considered with only two free parameters 
(these models treat only the in plane modes) are given by 
Kirkwood 6 , in which the elastic energy cost of stretching 
and bond bending are considered, and by Keating^, in a 
model that takes into account the symmetries of the lat- 
tice. These models can be extended to include the out of 
plane modes by adding a "dangling bond" term&£. The 
Kirkwood model is already a harmonic model, while the 
Keating model includes up to quartic order terms. The 
expansion up to second order of the Keating model gives 
an effective model very similar to the Kirkwood model, 
although it outperforms it slightly. For this reason in this 
paper we will take as a reference the extended Keating 
model — with the dangling bond term — as presented 
in Ref. |. 

In Ref. [l^ it was shown that for improving the ac- 
curacy of either the Kirkwood or Keating model, it is 
necessary to take into account the overlap of the 7r and 
a orbitals due to bending. In Ref. [l^ this was done by 



considering a full quantum mechanical model by use of a 
tight-binding formalism. On the other hand, in the work 
by Lenosky et al^, the authors propose a microscopic 
form for the energy of a graphene sheet that takes into 
account this overlap by considering the energy cost of 
having misaligned normal vectors to the graphene mem- 
brane. The original model was applied to treat the ener- 
getics of negatively curved graphene structures, denom- 
inated schwarzites and it was later extended to describe 
the elastic properties of nanotubes^. In this work we 
study this proposed free energy and calculate the resul- 
tant phonon dispersion, comparing with the available ex- 
perimental data. We will show that the terms involving 
the misalignment of the normals flatten the dispersion 
of the optical modes at the T point, a flattening that is 
observed experimentally — leaving aside, of course, the 
Kohn anomaly of the longitudinal optic (LO) model—!.. In 
particular the out-of-plane modes ZA (flexural acoustic) 
and ZO (flexural optic) are very well fitted. Graphene 
is the experimentally realized example of a polymerized 
membrane and these out-of-plane modes play a funda- 
mental role since it is the non-linear coupling between 
them and the in-plane modes which stabilizes the flat 
phased, and they are directly related to the rippling ob- 
served in the graphene sheets^. We will see however 
that, even though Lenosky's energy describes qualita- 
tively the phonon dispersion near the L point, it gives 
the wrong ordering in energy for the modes at the K 
point. We circumvent this by adding an extra term to 
Lenosky's original formulation, the bond bending term, 
which is present in both Kirkwood and Keating models 
but not in Lenosky's. 



II. THE MODEL 

In graphene the carbon atoms form a honeycomb lat- 
tice as shown in Fig. [TJ This configuration is due to the 
symmetry of the bonds between the carbon atoms: their 
s and p orbitals hybridize in the form sp 2 which results 
in cr bonds contained in the graphene plane, while the 
remaining p z orbital is perpendicular to the plane and 
forms the ir and ir* bands^. The proposed form for the 
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FIG. 1: Graphene lattice with the notation used throughout 
the text. 



< ij > denotes nearest neighbors. According to ReflTTI. 
the first two terms in expression |T]) correspond to bond 
stretching and angle bending respectively. We will how- 
ever refer to this second term as the "dangling bond" 
term since it is essentially the same term that was intro- 
duced as such in Ref. 0— . The last two terms take into 
account the energy cost of orbital overlapping due to rip- 
pling and are "non-local" in the sense that they involve 
more than nearest-neighbor terms as we will see below. 
The first of these two terms can be though as related to 
the 7r-7r orbital overlap, since it involves the scalar prod- 
uct between normal vectors at neighboring atoms, while 
the second is the projection of the normal onto the bonds 
and therefore related to the a-ir orbital overlap. 



energy given by Lenosky and collaborators in Ref. |ll| is 
given by the following expression: 
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where is the vector that points from atom i to atom 
j in the lattice, f ij — r^-r and ft; is the normal vector 

to the plane determined by the three nearest neighbors 
of atom i (we will call this loosely the "normal at atom 
i"). Vij = e.y for the undeformed lattice, where is 
a unit vector and we have absorbed the lattice constant 
a = 1.42 A in the definition of the elastic constants e„. 
The indices i, j run over the N atoms of the lattice and 



III. DYNAMICAL MATRIX 

For obtaining the phonon modes we derive the dynam- 
ical matrix determined by ([T]) via a harmonic expansion. 
If we allow for small displacements of the atoms with 
respect to their equilibrium positions in the lattice, we 
can write ry = + uj — uf , where uf is a small dis- 
placement of the i th atom in sublattice X, with X = A, 
B, and similarly for uj — . Up to second order in the 
displacements we get that expression (TTJ) can be written 
as U L = Wj° + Ul 1 + U e L 2 + , as we detail below. The 
stretching energy is given by the usual expression 
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The dangling bond term can be written in compact form 
by use of the dyadic notation as 
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We will denote from now on nearest neighbors with Greek 
indices. With this notation, indicates the displace- 

ment of the a th neighbor (a = 1, 2, 3) of atom n. The 
supra-index indicates that the displacement corresponds 
to an atom in sublattice X. The unit vectors i a and j Q 
are as indicated in Fig. [T]— and Id is the 3x3 identity 
matrix. 

The expansion for the terms that involve the normals 
to the graphene plane is more involved. Assuming a la- 
beling as in Fig. [T]we can write 



which is a vector of length equal to the area of a unit 
cell and it is normal to the n th atom of sublattice A, and 
analogous expressions hold for n^^. Using the explicit 
expressions for r$j in terms of the atoms' small displace- 
ments it can be shown that 

4 {B) = \^ + Y W >< « ) ( ; ) V3j 1 x uftf >) , 

where the sum is over cyclic permutations of the nearest 
neighbor index. Therefore the unit vector normal to atom 

n in sublattice X is given by nf — j^rj ■ It is straight- 
forward to retain the quadratic terms of expressions of 
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the type hf ■ hj , however for (fif • hf) we have 
to use a multivariable Taylor expansion. The expansion 

depends on nine variables since jn^l^ 1 - ' " x ,> x 1 
depends on the three 3D vectors u^m 
^ _ \„x\-i „ ro „„„ 



If we define the 



function f x (u^,^) = |n^ | 1 we can write 
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where is the m" 1 component of the small displace- 
ment of an atom in sublattice Y which is the neighbor \ 



of atom i in sublattice X. From the definition of the nor- 
mal vectors it is evident, as we mentioned previously, that 
terms that involve the product of normals at neighboring 
sites contain nevertheless displacements that go beyond 
nearest-neighbor interactions. Putting all together we 
obtain for the tt-tt and ir-a overlap energy cost the fol- 
lowing expressions, respectively: 
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(4) 
(5) 



From these expressions we can see that is related 
to acoustic modes while U^ 3 involves relative displace- 
ments of the A and B sublattices and then it is related 
to the optical modes. Both terms have components only 
in the out-of plane (z) direction and hence affect only the 
flexural modes. In-plane and out-of-plane modes are de- 
coupled in the harmonic approximation in the sense that 
the dynamical matrix can be block diagonalized. How- 
ever the modes are coupled through the elastic constant 
ei, since the dangling bond term has both in-plane 
and out-of-plane components. 



IV. RESULTS 

In Lenosky's original workii, the value of eo was taken 
to infinity and the remaining parameters were calculated 
by Local Density Approximation (LDA). The obtained 
values were e x = 0.96 eV, e 2 = 1.29 eV, e 3 = 0.05 eV. 
Here we take the parameter eo as adjustable to repro- 
duce correctly the value of the optical modes at the T 
point, from which we obtain eo = 36 eV. In Fig. [5] it is 
shown the phonon dispersion resulting from the Lenosky 
model plotted as a function of momentum with the above 
values for the e n parameters. As we anticipated, expres- 
sion (TT]) gives a very good description for the flexural 
modes near the T point, as it can be seen from Fig. [5] 
However the description is not so good for the in plane 
modes. The failing of the model for the LO mode at 
the r point and the transverse optic (TO) mode at the 
K point can be attributed to the Kohn anomalie o 20 ' 21 , 
that are not taken into account in |T]). However the most 
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FIG. 2: Solid line: phonon dispersion for the Lenosky model 
Vs. momentum with parameters as discussed in the text, 
dots: collection of experimental data taken from the review 
article Ref. @- The frequencies are in meV. 



important qualitative failure of the model is for the or- 
dering of the modes at the K point. With Lenosky's 
energy, the ZA/ZO modes at the K point have larger 
energy than the transverse acoustic (TA) mode, on the 
contrary of what is observed experimentally. A quick in- 
spection at the expressions of these modes in terms of 
the parameters e n reveals that this problem cannot be 
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FIG. 3: Solid line: best fit phonon dispersion for the hybrid 
model Vs. momentum with parameters as discussed in the 
text, dots: collection of experimental data taken from the 
review article Ref. 2 and Ref. 0. The frequencies are in meV. 



FIG. 4: Solid line: phonon dispersion for the hybrid model 
Vs. momentum resulting for a restricted fit as discussed in 
the text, dots: collection of experimental data taken from the 
review article Ref. 2 and Ref. [J The frequencies are in meV. 



solved by choosing a different set of parameters. The 
frequencies of the ZA/ZO and the TA modes at the K 
point are given by ujza{K) = \/a^f{\%€\ + 12e 2 ) JM 
and ujta(K) = \ j a\J\%t\jM respectively, being M the 
carbon mass, and where we have used the fact that 
eo has to be the largest energy scale in the problem. 
Since e„ > 0, from these expressions it is evident that 
lo Z a{K) > u TA (K). 

To overcome this issue we construct a "hybrid" model 
by adding to Lenosky's energy (TTJ) a bond-bending term 
as present in Keating and Kirkwood models: 

N , x x 2 

Ub = 22 zl 0(*ia(i)-ri ui i ) + -J , (6) 

i— 1 a, v~>ct ^ ' 

which represents the energy cost of changing the angle 
between bonds. This is roughly equivalent then to taking 
the model presented in Ref. and adding the it-it and tt-<t 
overlap terms, besides the small difference in the dangling 
bond term discussed previously. The harmonic expansion 
Ug of this term can be readily obtained and is given by 

U b = P 1 ' ( u ^(«) - u « ) + * 2 ' ( U U«) _ u n ) 

. (7) 

where the second sum is given over cyclic permutations 
of the nearest-neighbor indices. 

We present the results from our hybrid model in Fig. 
[3] for the high symmetry lines of the full Brioullin zone. 
We obtained a best fit for the experimental data in the 
T-K direction and utilized the resulting values for the 
elastic parameters e = 24.8 eV, t\ — 1.3 eV, ti = 0.2 
eV, £3 = 1.2 eV and (3 = 5.0 eV to construct the disper- 
sion for the rest of the Brioullin zone. From Fig. [3] it 
is seen that the fitting for the flexural modes is indeed 
exceptionally good. However problems still remain for 
the in-plane modes. In particular the degeneracy of the 



longitudinal acoustic (LA)/LO and TA/TO modes at the 
K point is not the right one. This problem is a conse- 
quence of the bond-bending tevmhlp and it is also present 
in both nearest-neighbor Keating and Kirkwood models, 
although it is not mentioned in the literature. The way 
out of this problem is to realize that these models allow 
for the right symmetry if a condition among the elastic 
parameters is fullfilled. We found this condition for our 
hybrid model to be eo > (7/3+ |ei) — which reduces 
to eo > 7/3 for the usual Keating model. Therefore, the 
fitting of the elastic parameters has to be constrained by 
this condition, instead of an unrestricted fit. If we apply 
this restriction to our fitting of the data, we obtain the 
phonon spectrum depicted in Fig. 21 From the figure 
it can be seen that the new fit has the right symmetries 
and it is overall a very good fit in the whole Brioullin 
zone, with the exception of the points where the Kohn 
anomaly should play a role: namely the overbending of 
the LO mode at the T point and the softening of the TO 
mode at the K point. From the values obtained for the 
elastic constants eo = 30 eV, ei = 1.3 eV, e 2 = 0.2 eV, 
£3 = 1.2 eV and (3 = 2.4 it can be inferred that the ir-a 
overlap dominates over the tt-tt orbital overla p. This is 
in agreement with the results presented in Ref. llOl within 
a tight binding study but differ with the values obtained 
for the Lenosky energy in Ref. 11 in which £3 is the 
smallest parameter. 

V. CONCLUSIONS 

In conclusion, we have analyzed the phonon disper- 
sion give n by the energy function |T]) first introduced in 
Ref. [ll| for graphene sheets. We have shown that this 
model gives good results for the flexural modes but fails 
to describe correctly in-plane modes beyond the T point. 
To overcome this we have constructed a five parameters 
hybrid model which adds to Lenosky's energy a com- 
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monly used bond-bending term. We have shown that 
a restricted fit of the elastic parameters reproduces cor- 
rectly all the relevant features of the phonon dispersion 
of graphene, with the exception of the Kohn anomalies. 
This restricted fit is essential to obtain the right symme- 
tries of the model, and our results with respect to this 
point also apply for existent and well established nearest- 
neighbor valence force models in the literature, namely, 
the Kirkwood 6 and Keating 7 models. According to our 
results, the effect of the change in overlap between the tt 
and a bonds is dominant over the effect due to the change 
in the it-it overlap. The agreement of our model with 
the experimental data is extremely good. In particular 
the flexural modes, which are of extreme importance for 



graphene, are excellently described. We have therefore 
obtained a minimal, few parameters model which can be 
used as a starting point in future analytical calculations. 
Results in this respect will be published elsewhere. 
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